% ***********************************************************************************************
% m_landowner.m 
% ***********************************************************************************************
clear all; 
close all;
format shortG;
% *************************************************************************
% ******************** LOAD FUNDAMENTAL PARAMETERS ************************
% *************************************************************************
load parameters.mat;                           
load calib_productivity_amenity.mat;     
load calib_setup.mat;                    
load commuting_cost.mat;  
load calib_auxiliary.mat;

load fundamentals_1950.mat;

load calib_qadjusted.mat; 
load calib_hscale.mat; 

Rdata=Rdata_all; Ldata=Ldata_all;  
totalpop=sum(Rdata,1);
N=size(ID,1);
NT=6;                       

% *************************************************************************
% *************************************************************************
% ****** Population in 1936 (pre-war) and determine distribution of *******
% ****** survived landowners based assuming uniform distribution of them **
% *************************************************************************
% *************************************************************************
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv'; 
CITYDATA=readtable(datafile);  
R36=CITYDATA.pop36; 
Pop36=sum(R36); 
ownership_rate=0.2;             
R36_owner=ownership_rate.*R36;
% survival rate data: From S22 Hiroshima shisei youran (April 1946)
survival_rate=zeros(N,1);
survival_rate=survival_rate+0.05.*(cbddist <= 1.0); 
survival_rate=survival_rate+0.14.*(1.0 < cbddist & cbddist <=1.5);
survival_rate=survival_rate+0.26.*(1.5 < cbddist & cbddist <=2.0);
survival_rate=survival_rate+0.80.*(2.0 < cbddist & cbddist <=2.5);
survival_rate=survival_rate+1.0.*(2.5 < cbddist);
% survival owners
R36_owner=survival_rate.*R36_owner; 
Pop36_owner=sum(R36_owner); 

% *************************************************************************
% ******************* COMPUTE COMMUTING PATTERN 1945 **********************
% *************************************************************************
Lcom=zeros(N,N,NT+1); 
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-5,'TolX',1e-5);
R45=Rdata(:,1); L45=Ldata(:,1); 
K45=Kappa_mat(:,:,1).^(-rho/sigma); 
Y0=ones(N,1);     
comsyst = @(y)fcomeval2(y,K45,R45,L45,rho,sigma);
[y_fsolve, yfval_fsolve]=fsolve(comsyst, Y0, options0); Y45=abs(y_fsolve); Y45=Y45./geomean(Y45);  
ytemp=K45.*(repmat(Y45,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
Lcom45=ytemp.*repmat(R45',N,1); 
Lcom(:,:,1)=Lcom45;

% *************************************************************************
% ********************* EXOGENOUS FUNDAMENTALS (50-75) ********************
% ********************* OPTION VALUES OF FUNDAMENTALS *********************
% *************************************************************************
a_mat=[a50, aa];   a_mat=a_mat./repmat(geomean(a_mat),N,1);
b_mat=[b50, bb];   b_mat=b_mat./repmat(geomean(b_mat),N,1);
Xmat_noagg=zeros(N,NT); Ymat_noagg=zeros(N,NT); 
Xmat_noagg(:,NT)=b_mat(:,NT); Ymat_noagg(:,NT)=a_mat(:,NT);
for t=1:NT-1
    tt = NT-t; 
    theta=thetavec(tt);
    Xtemp=Xmat_noagg(:,tt+1).^(rho*(1-theta)); Xtemp=b_mat(:,tt).*Xtemp;
    Ytemp=Ymat_noagg(:,tt+1).^(rho*(1-theta)); Ytemp=a_mat(:,tt).*Ytemp;
    Xmat_noagg(:,tt)=Xtemp;
    Ymat_noagg(:,tt)=Ytemp;
end
% *************************************************************************
% *************************************************************************
% ***************** COMPUTE COUNTERFACTUAL FROM 1950 **********************
% *************************************************************************
% *************************************************************************
Q_i=ones(N,NT); 
L_i=Ldata(:,2:end); R_i=Rdata(:,2:end);  

Q_n=zeros(N,NT); L_n=zeros(N,NT); R_n=zeros(N,NT); H_n=zeros(N,NT+1); H_n(:,1)=hsupply45; 
H_demand=zeros(N,NT); 

maxit=1e+6; tol=1e-6; update=0.05; it=0; 

while it<maxit  
   
LQmat=zeros(N,NT); 
LQmat(:,NT)=Q_i(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    LQmat(:,tt)=Q_i(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta)));
end 

% main computation from 1950 to 1975 
for t=1:NT

    if t==1; theta=theta50; Rpre_i=R45; Lpre_i=L45; 
    else     theta=thetavec(t-1); Rpre_i=R_n(:,t-1); Lpre_i=L_n(:,t-1); end 
    
    % conditional probabilities (lambda_L: residential given workplace; lambda_R: workplace given residenace) 
    lambda_L=Kmat(:,:,t).*((repmat(LQmat(:,t)',N,1).^(-mu)).*(repmat(Xmat_noagg(:,t)',N,1))).^(rho/sigma); 
    lambda_L=lambda_L./repmat(sum(lambda_L,2),1,N); 
    lambda_R=Kmat(:,:,t).*((repmat(LQmat(:,t),1,N).^(-gammaH/gammaL)).*(repmat(Ymat_noagg(:,t),1,N).^(1/gammaL))).^(rho/sigma);
    lambda_R=lambda_R./repmat(sum(lambda_R,1),N,1);

    % initial guess of population/employment 
    % survived owner population + population share *  Non-ownerpopulation for 1950
    if t==1 
        Lpost_i=R36_owner+(Ldata(:,1)./totalpop(1)).*(totalpop(2)-Pop36_owner); 
        Rpost_i=R36_owner+(Rdata(:,1)./totalpop(1)).*(totalpop(2)-Pop36_owner); 
    else
        Rpost_i=R_i(:,t); Lpost_i=L_i(:,t);  
    end 
    
    it1=0;

    while it1<maxit 

    if t==1
        dR=Rpost_i-(1-theta).*R45-R36_owner; dR=abs(dR); 
        dL=Lpost_i-(1-theta).*L45-R36_owner; dL=abs(dL); 
        Rpost_n=R36_owner+(1-theta).*Rpre_i+lambda_L'*dL;  
        Lpost_n=R36_owner+(1-theta).*Lpre_i+lambda_R*dR;
    else 
        Rpost_n=(1-theta).*Rpre_i+lambda_L'*(Lpost_i-(1-theta).*Lpre_i);
        Lpost_n=(1-theta).*Lpre_i+lambda_R*(Rpost_i-(1-theta).*Rpre_i);  
    end 

    if max(abs(Rpost_n-Rpost_i))<tol && max(abs(Lpost_n-Lpost_i))<tol; break;
    else
       Rpost_i=update.*Rpost_n+(1-update).*Rpost_i; 
       Lpost_i=update.*Lpost_n+(1-update).*Lpost_i;
       it1=it1+1;  
    end
    end
    R_n(:,t)=Rpost_n; L_n(:,t)=Lpost_n; 

    % commuting patterns (column:=workplace; row:=residential place)
    if t==1
        Lcompost=diag(R36_owner)+(1-theta).*Lcom(:,:,t)+lambda_R.*repmat(dR',N,1); % number of commuters
    else
        Lcompost=(1-theta).*Lcom(:,:,t)+lambda_R.*repmat((Rpost_n-(1-theta).*Rpre_i)',N,1); % number of commuters
    end
    Lcom(:,:,t+1)=Lcompost;
    
    % wage (workplace) and average income (residential place) 
    w_n=(Q_i(:,t).^(-gammaH/gammaL)).*(a_mat(:,t).^(1/gammaL));
    v_n=(Lcompost./repmat(Rpost_n',N,1))'*w_n; 

    % aggregate demand and supply (stock) of floor space
    h_n=(1-psi).*H_n(:,t)+hbar(t)*(Q_i(:,t).^eta).*Area;
    H_n(:,t+1)=h_n; 
    floor_demand=mu.*v_n.*Rpost_n + (gammaH/gammaL).*w_n.*Lpost_n;
    H_demand(:,t)=floor_demand; 
    
    % floor space market clearing condition
    Q_n(:,t)=floor_demand./h_n ;
end 

dd_Q=Q_n-Q_i; 

if max(max(abs(dd_Q))) < tol 
        disp(['Converged floor prices']); break;
else
   Q_i=update.*Q_n+(1-update).*Q_i; it=it+1;  
   if rem(it,10)==1
      disp(['Largest Error is:',' ',num2str(it), ' ',num2str(max(abs(dd_Q)))]);
   end
end

end

% check floor space per area, as well as the ratio of existing stock and
% total floor space
max(H_n./Area)
median(H_n./Area)

% check floor market clearing 
if max(abs(Q_n.*H_n(:,2:end)-H_demand))<1e-6; disp(['PASS floor space market clearing conditions in 1950-75']);
else; disp(['WARNING']); end 

% *************************************************************************
% *************************** SAVE RESULTS ********************************
% *************************************************************************
R_noagg=R_n; L_noagg=L_n; Q_noagg=Q_n; H_noagg=H_n; Lcom_noagg=Lcom; 

result=array2table([R_noagg,Rdata,L_noagg,Ldata, Q_noagg, Qadj, hstock, H_noagg, cbddist, Area, ID, lat, lon],'VariableNames',{'pop50_noagg',...
    'pop55_noagg','pop60_noagg','pop65_noagg','pop70_noagg','pop75_noagg',...
    'pop45','pop50','pop55','pop60','pop65','pop70','pop75',...
    'emp50_noagg','emp55_noagg','emp60_noagg','emp65_noagg','emp70_noagg','emp75_noagg',...
    'emp45','emp50','emp55','emp60','emp65','emp70','emp75',...
    'q50_noagg','q55_noagg','q60_noagg','q65_noagg','q70_noagg','q75_noagg',... 
    'q50','q55','q60','q65','q70','q75',... 
    'h50','h55','h60','h65','h70','h75',... 
    'h45','h50_noagg','h55_noagg','h60_noagg','h65_noagg','h70_noagg','h75_noagg',...
    'cbddist', 'area', 'geocode', 'lat', 'lon'}); 
writetable(result,'../../output/quant/output_no_agglomeration_landowner.csv');
% *************************************************************************
% ************************** FUNCTIONS ************************************ 
% *************************************************************************
function[comeval] = fcomeval2(Y,K,DR,DL,rho,sigma) 

N = size(DR,1); 
Y=abs(Y);  
Y=Y./geomean(Y);
ytemp=K.*(repmat(Y,1,N).^(rho/sigma)); 
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 
comeval = DL - sum(ytemp.*repmat(DR',N,1),2); 
comeval = abs(comeval); 

end

